/**
This file is part of Gel2DE.

    Gel2DE is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Gel2DE is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with Gel2DE.  If not, see <http://www.gnu.org/licenses/>.

	Copyright (C) 2012, University of Bergen
*/

#include "GelCorrelationAnalysis.h"
#include "SpearmanRankOrderCorrCoeff.h"
#include "tinyxml.h"
#include <wx/wx.h>
#include <wx/progdlg.h>

#include "vtkColorTransferFunction.h"

#include <stdlib.h>
#include <time.h>
#include <sstream>
#include <iostream>
#include <fstream>

GelCorrelationAnalysis::GelCorrelationAnalysis(vtkFileOutputWindow *aLogFile)
{
	m_ctf = NULL;
	m_gelVector = new std::vector<GelImage*>();
	m_gelImageTransformationMap = new std::map<const std::string, GelImageTransformation*>();
	m_logFile = aLogFile;
	srand(time(NULL));
}

GelCorrelationAnalysis::~GelCorrelationAnalysis()
{
}


void GelCorrelationAnalysis::Clear()
{
	std::vector<GelImage*>::iterator iter = m_gelVector->begin();
	for(iter ; iter != m_gelVector->end(); ++iter)
	{
		delete (*iter);
	}

	m_gelVector->clear();
	m_gelImageTransformationMap->clear();
}

void GelCorrelationAnalysis::SetBlurLevel(double aSigma)
{
	std::vector<GelImage*>::iterator i = m_gelVector->begin();
	for(i; i != m_gelVector->end(); ++i)
	{
		(*i)->SetBlurLevel(aSigma);
	}
}
/*
  Normalization mode: 0 - None, 1 - Mean, 2 - Median
*/
void GelCorrelationAnalysis::SetNormalizationMode(int aNormalizationMode)
{
	GelImage* frontImage = m_gelVector->front();
	double globalScale = 1.;
	if (aNormalizationMode == GelImage::MEAN_NORM) globalScale = frontImage->GetAvgIntensity();
	else if(aNormalizationMode == GelImage::MEDIAN_NORM) globalScale = frontImage->GetMedian();

	std::vector<GelImage*>::iterator i = m_gelVector->begin();
	for(i; i != m_gelVector->end(); ++i)
	{
		(*i)->SetNormalizationMode(aNormalizationMode, globalScale);
	}
}


std::vector<GelImage*>* GelCorrelationAnalysis::GetGelImages()
{
	return m_gelVector;
}

GelImage* GelCorrelationAnalysis::GetImage(int aImageNumber)
{
	if(aImageNumber < m_gelVector->size())
	{
		return (*m_gelVector)[aImageNumber];
	}
	else
	{
		return NULL;
	}
}

void GelCorrelationAnalysis::LoadGelImageTransformationMap(std::string aFileName)
{
	TiXmlNode* stackNode;

	m_gelImageTransformationMap->clear();

	TiXmlDocument doc(aFileName.c_str());
	if(!doc.LoadFile())
	{
		std::string errorMsg("WARNING: Registration file " + aFileName + " not found!");
		m_logFile->DisplayText(errorMsg.c_str());
	}
	else
	{
		std::string infoMsg("Reading " + aFileName + "...");
		m_logFile->DisplayText(infoMsg.c_str());
	}
	for(stackNode = doc.FirstChild(); stackNode != 0; stackNode = stackNode->NextSibling())
	{
		std::string name = stackNode->Value();
		if(stackNode->Type() == TiXmlNode::TINYXML_ELEMENT && name == "Stack")
		{
			TiXmlNode* gelImageNode;
			for(gelImageNode = stackNode->FirstChild(); gelImageNode != 0; gelImageNode = gelImageNode->NextSibling())
			{
				name = gelImageNode->Value();
				if(gelImageNode->Type() == TiXmlNode::TINYXML_ELEMENT && name == "GelImage")
				{
					std::string aFileName;
					TiXmlNode* childNode;
					for(childNode = gelImageNode->FirstChild(); childNode != 0; childNode = childNode->NextSibling())
					{
						name = childNode->Value();
						if(childNode->Type() == TiXmlNode::TINYXML_ELEMENT && name == "FileName")
						{
							aFileName = ((TiXmlElement*)childNode)->GetText();
						}
					}

					GelImageTransformation* transf = new GelImageTransformation((TiXmlElement*)gelImageNode);

					m_gelImageTransformationMap->insert(std::pair<const std::string, GelImageTransformation*>(aFileName, transf));
				}
			}
		}
	}
}

void GelCorrelationAnalysis::LoadImageAndClinicalParameters(std::string aFileName, std::vector<std::string*>* aExcludedList)
{
	TiXmlNode* stackNode;

	TiXmlDocument doc(aFileName.c_str());
	if(!doc.LoadFile())
	{
		std::string msg("ERROR: " + aFileName + " could not be parsed!");
		m_logFile->DisplayText(msg.c_str());
	}
	else
	{
		std::string msg("Reading " + aFileName + "...");
		m_logFile->DisplayText(msg.c_str());
	}
	for(stackNode = doc.FirstChild(); stackNode != 0; stackNode = stackNode->NextSibling())
	{
		std::string name = stackNode->Value();
		if(stackNode->Type() == TiXmlNode::TINYXML_ELEMENT && name == "Stack")
		{
			TiXmlNode* gelImageNode;
			int numberofGels = 0;
			for(gelImageNode = stackNode->FirstChild(); gelImageNode != 0; gelImageNode = gelImageNode->NextSibling())
			{
				numberofGels++;
			}

			wxProgressDialog *progressDlg = new wxProgressDialog(wxT(""), wxT("Loading project, please wait..."), numberofGels);
			int numberOfLoadedGels = 0;

			for(gelImageNode = stackNode->FirstChild(); gelImageNode != 0; gelImageNode = gelImageNode->NextSibling())
			{
				numberOfLoadedGels++;
				progressDlg->Update(numberOfLoadedGels);
				name = gelImageNode->Value();
				if(gelImageNode->Type() == TiXmlNode::TINYXML_ELEMENT && name == "GelImage")
				{
					std::string gelImageFileName;
					TiXmlNode* childNode;
					for(childNode = gelImageNode->FirstChild(); childNode != 0; childNode = childNode->NextSibling())
					{
						//Now looping over the GelImage tags. Find the filename to look up the transformation
						name = childNode->Value();
						if(childNode->Type() == TiXmlNode::TINYXML_ELEMENT && name == "FileName")
						{
							gelImageFileName = ((TiXmlElement*)childNode)->GetText();
						}
					}
					std::map<const std::string, GelImageTransformation*>::iterator iter = m_gelImageTransformationMap->find(gelImageFileName);
					GelImage *gelImage;
					if(iter != m_gelImageTransformationMap->end())
					{
						gelImage = new GelImage((TiXmlElement*)gelImageNode, iter->second, m_logFile);
					}
					else
					{
						gelImage = new GelImage((TiXmlElement*)gelImageNode, NULL, m_logFile);
					}
					//Is this enabled:
					std::vector<std::string*>::iterator iexcluded = aExcludedList->begin();
					for(iexcluded; iexcluded != aExcludedList->end(); ++iexcluded)
					{
						if((*(*iexcluded)) == gelImage->GetFileName())
						{
							gelImage->SetEnabled(FALSE);
							break;
						}
					}
					
					m_gelVector->push_back(gelImage);

				}
			}
			delete progressDlg;
		}
	}
}

void GelCorrelationAnalysis::ExportData(GelCorrelationAnalysis::RGBImageType::RegionType aRoi, std::string aFileName)
{

	ofstream outputfile;
	outputfile.open(aFileName.c_str());

	std::vector<GelImage*>::iterator i = m_gelVector->begin();

	int numberOfGels = m_gelVector->size();
	wxProgressDialog *progressDlg = new wxProgressDialog(wxT(""), wxT("Exporting data, please wait..."), numberOfGels);
	int numberOfLoadedGels = 0;

	GelImage::ImageType::RegionType gelRoi;
	//DataRoiToGelRoi(&aRoi, &gelRoi);

	for(i; i != m_gelVector->end(); ++i)
	{
		numberOfLoadedGels++;
		progressDlg->Update(numberOfLoadedGels);
		GelImage* image = (*i);
		image->SetRoi(aRoi);
		image->MoveToStart();
		outputfile << image->GetFileName() << std::endl;
		image->Update();
		while (!image->IsAtEnd())
		{
			outputfile << image->GetCurrentPixelRaw() << ";" ;
			image->MoveToNextPixel();
		}
		outputfile << std::endl;
		std::set<const std::string> parameterNames = image->GetParameterTypes();
		std::set<const std::string>::iterator parIter = parameterNames.begin();
		for( parIter; parIter != parameterNames.end(); ++parIter)
		{
			outputfile << (*parIter) << ";" << image->GetParameterValue((*parIter)) << std::endl;
		}
		outputfile << std::endl;
	}
	outputfile.close();
	delete progressDlg;
}

void GelCorrelationAnalysis::SetOutputDirectory(gel2de::File* aFile)
{
	m_outputDirectory = aFile;
}


void GelCorrelationAnalysis::LoadConfiguration(ProjectConfiguration* aConfiguration)
{

	LoadGelImageTransformationMap(aConfiguration->getRegistrationFile()->name());

	LoadImageAndClinicalParameters(aConfiguration->getParameterFile()->name(), aConfiguration->getExcludedFilesStoredInProject());
//todo: dangerous giving in the configuration object directly?
	SetOutputDirectory(aConfiguration->getOutputDataDir());
}

std::set<const std::string> GelCorrelationAnalysis::getAllParameterTypes()
{
	std::set<const std::string> allTypes;
	std::vector<GelImage*>::iterator gelImageIter = m_gelVector->begin();
	for(gelImageIter; gelImageIter != m_gelVector->end(); ++gelImageIter)
	{
		std::set<const std::string> types = (*gelImageIter)->GetParameterTypes();

		std::set<const std::string> temp_set;
		std::set_union(
			allTypes.begin(), allTypes.end(), 
			types.begin(), types.end(), 
			std::inserter(temp_set, temp_set.begin())
			);
		allTypes.swap(temp_set);
	}

	return allTypes;
}

void GelCorrelationAnalysis::ToXML(const char* aFileName)
{
	TiXmlDocument doc;
	TiXmlDeclaration* decl = new TiXmlDeclaration("1.0", "", "");
	TiXmlElement* element = new TiXmlElement("Stack");
	std::vector<GelImage*>::iterator iter = m_gelVector->begin();
	for(iter; iter != m_gelVector->end(); iter++)
	{
		TiXmlElement* gelImageElement = (*iter)->ToXML();
		element->LinkEndChild(gelImageElement);
	}
	doc.LinkEndChild(decl);
	doc.LinkEndChild(element);
	doc.SaveFile(aFileName);
}

GelCorrelationAnalysis::RGBPixelType* GelCorrelationAnalysis::ToRGB(double aValue)
{
  double r = m_ctf->GetRedValue(aValue);
  double g = m_ctf->GetGreenValue(aValue);
  double b = m_ctf->GetBlueValue(aValue);

  RGBComponentType rgb[3] = {r*255, g*255, b*255};
  return new RGBPixelType(rgb);
}


GelCorrelationAnalysis::RGBImageType::Pointer GelCorrelationAnalysis::createRGBImage(
				GelImage::ImageType::RegionType aSourceRegion,
				const GelImage::ImageType::SpacingType& aSpacing,
				const GelImage::ImageType::PointType& aInputOrigin)

{

	RGBImageType::RegionType outputRegion;
	RGBImageType::Pointer outImg = RGBImageType::New();

	outputRegion.SetSize(aSourceRegion.GetSize());
	outputRegion.SetIndex(aSourceRegion.GetIndex());

	outImg->SetRegions(outputRegion);

	outImg->SetSpacing(aSpacing);
	outImg->SetOrigin(aInputOrigin);
	outImg->Allocate();

	return outImg;
}

GelCorrelationAnalysis::DataImageType::Pointer GelCorrelationAnalysis::createDataImage(
	GelImage::ImageType::RegionType aSourceRegion,
	const GelImage::ImageType::SpacingType& aSpacing,
	const GelImage::ImageType::PointType& aInputOrigin)

{
	DataImageType::RegionType outputRegion;
	DataImageType::Pointer outImg = DataImageType::New();

	outputRegion.SetSize(aSourceRegion.GetSize());
	outputRegion.SetIndex(aSourceRegion.GetIndex());

	outImg->SetRegions(outputRegion);

	outImg->SetSpacing(aSpacing);
	outImg->SetOrigin(aInputOrigin);
	outImg->Allocate();

	return outImg;
}

int GelCorrelationAnalysis::GetNumberOfActiveGels(std::string aParameter)
{
	std::vector<GelImage*>::iterator igel;

	int N = 0;
	igel = m_gelVector->begin();
	for(igel ; igel != m_gelVector->end(); ++igel)
	{
		if((*igel)->isEnabled() && (*igel)->HasParameterValue(aParameter))
		{
			N++;
		}
	}

	return N;
}

double GelCorrelationAnalysis::GetAverage(std::string aParameter)
{
	std::vector<GelImage*>::iterator igel;

	double sum = 0;
	int N = 0;
	igel = m_gelVector->begin();
	for(igel ; igel != m_gelVector->end(); ++igel)
	{
		if((*igel)->isEnabled() && (*igel)->HasParameterValue(aParameter))
		{
			sum += (*igel)->GetParameterValue(aParameter);
			N++;
		}
	}

	return sum/N;
}

void GelCorrelationAnalysis::DataRoiToGelRoi(DataImageType::RegionType *aDataRoi, GelImage::ImageType::RegionType *aGelRoi)
{
	DataImageType::RegionType::IndexType index = aDataRoi->GetIndex();
	DataImageType::RegionType::SizeType size = aDataRoi->GetSize();
	GelImage::ImageType::RegionType::IndexType rawIndex;
	GelImage::ImageType::RegionType::SizeType rawSize;

	rawIndex[0] = index[0];
	rawIndex[1] = index[1];
	rawSize[0] = size[0];
	rawSize[1] = size[1];

	aGelRoi->SetIndex(rawIndex);
	aGelRoi->SetSize(rawSize);
}

void GelCorrelationAnalysis::RunPermutationTest(DataImageType::RegionType aRoi, std::set<const std::string> aListOfParameters, bool aWriteToFile)
{
	//we should now make iterators for each image
	m_logFile->DisplayText("Setting up ROI");

	DataImageIteratorType ptIter(m_outImgPT, aRoi);
	ptIter.GoToBegin();

	GelImage::ImageType::RegionType roi;
	DataRoiToGelRoi(&aRoi, &roi);
	m_logFile->DisplayText("Set up ROI");

	//Now we can create iterators.
	std::vector<GelImage*>::const_iterator igel = m_gelVector->begin();

	int npixels = roi.GetNumberOfPixels();

	wxProgressDialog *progressDlg = new wxProgressDialog(wxT(""), wxT("Please be patient while performing the permutation test."),npixels);

	char out[1024];
	sprintf(out, "n pixels %d", npixels);
	m_logFile->DisplayText(out);
	for(igel = m_gelVector->begin(); igel != m_gelVector->end(); ++igel)
	{
		(*igel)->UseData();
		(*igel)->Update();
		(*igel)->SetRoi(roi);
		(*igel)->MoveToStart();
	}
	int ngels = m_gelVector->size();
	std::vector<double> intensities(ngels);
	std::vector<double> values(ngels);
	const int npermutations = 5000;
	std::vector<double> permutedValues(npermutations);
	m_logFile->DisplayText("Starting loop...");


	bool endOfImage = false;
	int currentpixel = 0;
	while(!endOfImage)
	{
		intensities.clear();
		values.clear();
		progressDlg->Update(currentpixel);
		currentpixel++;

		int i = 0;
		for(igel = m_gelVector->begin() ; igel != m_gelVector->end(); ++igel)
		{	
			if((*igel)->isEnabled() && (*igel)->HasParameterValue(*(aListOfParameters.begin())))
			{
				double I = ((double)(*igel)->GetCurrentPixelRaw());///(*igel)->GetAvgIntensity();
				double T = (*igel)->GetParameterValue((*(aListOfParameters.begin())));

				intensities.push_back(I);
				values.push_back(T);

				i++;
				(*igel)->MoveToNextPixel();
				if((*igel)->IsAtEnd())
					endOfImage = TRUE;
			}
		}
		//First, recalculate the unpermuted value:
		SpearmanRankOrderCorrCoeff *calc = new SpearmanRankOrderCorrCoeff(&intensities, &values);
		std::pair<double, double>* result = calc->calculateCorrCoeff(false);
		delete calc;
		double r_unpermuted = fabs(result->first);

		std::string fileName;
		std::stringstream sstream;
		DataImageType::IndexType currentIndex = ptIter.GetIndex();
		sstream << "column_" << currentIndex[0] << "_" << currentIndex[1] << ".txt";
		
		gel2de::File absFile((*m_previousAnalysisOutputDir), sstream.str());

		ofstream outfile;
		if (aWriteToFile)
		{
			outfile.open(absFile.path().c_str());

			outfile << "Stack column:" << std::endl;
			for (int i = 0; i < intensities.size(); i++)
			{
				outfile << intensities[i] << ";" << values[i] << std::endl;
			}
			outfile << "Measured corr: " << std::endl;
			outfile << result->first << std::endl;
			outfile << "Permuted sets:" << std::endl;
		}

		SpearmanRankOrderCorrCoeff *permutedCalc = new SpearmanRankOrderCorrCoeff(ngels);
		for(int i = 0; i < npermutations; i++)
		{
			permutedCalc->Clear();
			Permute(&values);
			permutedCalc->AddData(&intensities, &values);
			std::pair<double, double>* pair = permutedCalc->calculateCorrCoeff(false);
			permutedValues[i] = fabs(pair->first);
			if (aWriteToFile) outfile << pair->first << ";";
			delete pair;
		}
		if (aWriteToFile) outfile << std::endl;
		delete permutedCalc;

		int nsmaller = 0;
		int nlarger = 0;
		for(int i = 0; i < npermutations; i++)
		{
			if (permutedValues[i] < r_unpermuted)
			{
				nsmaller++;
			}
			else
			{
				nlarger++;
			}
		}

		char out[1024];
		double permutedP = (double)nlarger/(nlarger + nsmaller);
		ptIter.Set(permutedP);
		++ptIter;
		if (aWriteToFile)
		{
			sprintf(out, "R: %f, N larger %d smaller %d, p-value %f", result->first, nlarger, nsmaller, permutedP);
			m_logFile->DisplayText(out);

			outfile << "Permutation p-value:" << std::endl;
			outfile << permutedP << std::endl;
			outfile.close();
		}
	}

	for(igel = m_gelVector->begin(); igel != m_gelVector->end(); ++igel)
	{
		(*igel)->SetupIterator();
	}

	GelImage::ImageType::Pointer sourceImg = m_gelVector->front()->GetImage();
	GelImage::ImageType::RegionType region = sourceImg->GetRequestedRegion();
	const GelImage::ImageType::SpacingType& spacing = sourceImg->GetSpacing();
	const GelImage::ImageType::PointType& inputOrigin = sourceImg->GetOrigin();
	DataImageType::RegionType outputRegion = m_outImgC->GetRequestedRegion();
	DataImageIteratorType outIterPT(m_outImgPT, outputRegion);

	RGBImageType::Pointer imgRGBPT = FillRGBImage(region, spacing, inputOrigin,outputRegion, outIterPT);
	WriteResultFile("PT", m_idstring, m_previousAnalysisOutputDir, imgRGBPT);

	delete progressDlg;

}

void GelCorrelationAnalysis::RunAnalysis(std::set<const std::string> aListOfParameters)
{

	typedef itk::ImageRegionIterator<GelImage::ImageType> IteratorType;

	std::cout << " - Calculating correlation" << std::endl;

	GelImage::ImageType::Pointer sourceImg = m_gelVector->front()->GetImage();
	GelImage::ImageType::RegionType region = sourceImg->GetRequestedRegion();
	const GelImage::ImageType::SpacingType& spacing = sourceImg->GetSpacing();
	const GelImage::ImageType::PointType& inputOrigin = sourceImg->GetOrigin();

	int maxPixelCOunt = region.GetNumberOfPixels();


	int ngels = m_gelVector->size();

	//Should be moved to GUI code.
	wxProgressDialog *progressDlgPre = new wxProgressDialog(wxT(""), wxT("Preparing for calculation..."),ngels);

	std::vector<GelImage*>::iterator igel;

	int gelnumber = 0;
	for(igel = m_gelVector->begin(); igel != m_gelVector->end(); ++igel)
	{
		progressDlgPre->Update(gelnumber);
		gelnumber++;
		(*igel)->UseData();
		(*igel)->Update();
		(*igel)->SetupIterator();
	}

	m_outImgC   = createDataImage(region, spacing, inputOrigin);
	m_outImgS   = createDataImage(region, spacing, inputOrigin);
	m_outImgP   = createDataImage(region, spacing, inputOrigin);
	m_outImgPT  = createDataImage(region, spacing, inputOrigin);
	m_outImgD   = createDataImage(region, spacing, inputOrigin);
	m_outImgCS  = createDataImage(region, spacing, inputOrigin);
	m_outImgCSD = createDataImage(region, spacing, inputOrigin);

	DataImageType::RegionType outputRegion = m_outImgC->GetRequestedRegion();

	DataImageIteratorType outIterC(m_outImgC, outputRegion);
	DataImageIteratorType outIterS(m_outImgS, outputRegion);
	DataImageIteratorType outIterP(m_outImgP, outputRegion);
	DataImageIteratorType outIterPT(m_outImgPT, outputRegion);
	DataImageIteratorType outIterD(m_outImgD, outputRegion);
	DataImageIteratorType outIterCS(m_outImgCS, outputRegion);
	DataImageIteratorType outIterCSD(m_outImgCSD, outputRegion);

	outIterC.GoToBegin();
	outIterS.GoToBegin();
	outIterP.GoToBegin();
	outIterPT.GoToBegin();
	outIterD.GoToBegin();
	outIterCS.GoToBegin();
	outIterCSD.GoToBegin();

	delete progressDlgPre;
	bool endOfImage = FALSE;

	//Should be moved to GUI code.
	wxProgressDialog *progressDlg = new wxProgressDialog(wxT(""), wxT("Please be patient while calculating correlation!"),maxPixelCOunt);

	std::cout << " - Calculating correlations... " << std::endl;
	int pixelcount = 0;

	double maxD = -1;

	SpearmanRankOrderCorrCoeff* calculator = new SpearmanRankOrderCorrCoeff(ngels);

	while(!endOfImage)
	{
		calculator->Clear();
		if(pixelcount % 10000 == 0)
		{
			progressDlg->Update(pixelcount);
			std::cout << "Pixel " << pixelcount << std::endl;
		}
		pixelcount++;

		//Iterate over all gels
		igel = m_gelVector->begin();

		for(igel ; igel != m_gelVector->end(); ++igel)
		{
			if((*igel)->isEnabled() && (*igel)->HasParameterValue(*(aListOfParameters.begin())))
			{
				//Normalize pixel intensity with the average intensity in this image.
				double I = ((double)(*igel)->GetCurrentPixelRaw());///(*igel)->GetAvgIntensity();

				//Get the value of the clinical parameter (only first in list for now)
				double T = (*igel)->GetParameterValue((*(aListOfParameters.begin())));

				calculator->AddData(I, T);
				//Move the iterator of current gel to next pixel so that it is
				//positioned correctly for next pixel calculation
				(*igel)->MoveToNextPixel();

				//calculator->AddData((double)I, T);
				if((*igel)->IsAtEnd())
					endOfImage = TRUE;
			}
		}
		
		std::pair<double, double>* result = calculator->calculateCorrCoeff();

		double rs = result->first;

		DataPixelType pixelValC = rs;
		outIterC.Set(pixelValC);

		double S = 1 - result->second;
		DataPixelType pixelValS = S;
		outIterS.Set(pixelValS);

		double P = result->second;
		DataPixelType pixelValP = P;
		outIterP.Set(pixelValP);
		
		outIterPT.Set(-1);

		double stddev = calculator->StdDevI(); 
		DataPixelType pixelValD = stddev;
		outIterD.Set(pixelValD);

		if(pixelValD > maxD)
		{
			maxD = pixelValD;
		}

		++outIterC;
		++outIterS;
		++outIterP;
		++outIterPT;
		++outIterD;

		delete result;
	}

	//Move all iterators back to start for new calculation
	igel = m_gelVector->begin();
	for(igel ; igel != m_gelVector->end(); igel++)
	{
		(*igel)->MoveToStart();
	}

	//Find maximum stddev:
	NormalizeImage(m_outImgD, outputRegion, maxD);

	//Now that the stddev is normalized we can create combination images
	outIterC.GoToBegin();
	outIterS.GoToBegin();
	outIterD.GoToBegin();

	while(!outIterC.IsAtEnd())
	{
		outIterCS.Set(outIterC.Get()*outIterS.Get());
		outIterCSD.Set(outIterC.Get()*outIterS.Get()*outIterD.Get());

		++outIterC;
		++outIterS;
		++outIterD;
		++outIterCS;
		++outIterCSD;
	}

	RGBImageType::Pointer outRGBImgC = FillRGBImage(region, spacing, inputOrigin, outputRegion, outIterC);
	RGBImageType::Pointer outRGBImgCS = FillRGBImage(region, spacing, inputOrigin, outputRegion, outIterCS);
	RGBImageType::Pointer outRGBImgCSD = FillRGBImage(region, spacing, inputOrigin, outputRegion, outIterCSD);
	RGBImageType::Pointer outRGBImgS = FillRGBImage(region, spacing, inputOrigin, outputRegion, outIterS);
	RGBImageType::Pointer outRGBImgD = FillRGBImage(region, spacing, inputOrigin, outputRegion, outIterD);
	RGBImageType::Pointer outRGBImgP = FillRGBImage(region, spacing, inputOrigin, outputRegion, outIterP);
	RGBImageType::Pointer outRGBImgPT = FillRGBImage(region, spacing, inputOrigin, outputRegion, outIterPT);

	std::string base("corr");
	char dateStr[9];
	char timeStr[9];
	_strdate_s(dateStr);
	_strtime_s(timeStr);

	std::string::size_type k = 0;

	std::string date(dateStr);
	while( (k = date.find('/',k)) != date.npos) 
	{
		date.erase(k, 1);
	}

	k = 0;
	std::string time(timeStr);
	while( (k=time.find(':',k)) != time.npos) 
	{
		time.erase(k, 1);
	}

	double blurLevel = m_gelVector->front()->GetBlurLevel();
	char blurText[64];
	sprintf(blurText, "b%d", int(blurLevel*10));

	int normType = m_gelVector->front()->GetNormalizationMode();
	char normText[64];
	if (normType == GelImage::NO_NORM) sprintf(normText, "no");
	else if (normType == GelImage::MEAN_NORM) sprintf(normText, "mean");
	else if (normType == GelImage::MEDIAN_NORM) sprintf(normText, "med");
	
	std::string blurStr(blurText);
	std::string normStr(normText);

	std::string idstring(date + "-" + time + "-" + normStr + "-" + blurStr + "-" + *(aListOfParameters.begin()));

	gel2de::File* path = new gel2de::File((*m_outputDirectory), std::string("corr-") + idstring);

	_mkdir(path->path().c_str());

	WriteResultFile("C", idstring, path, outRGBImgC);
	WriteResultFile("CS", idstring, path, outRGBImgCS);
	WriteResultFile("CSD", idstring, path, outRGBImgCSD);
	WriteResultFile("S", idstring, path, outRGBImgS);
	WriteResultFile("D", idstring, path, outRGBImgD);
	WriteResultFile("P", idstring, path, outRGBImgP);
	WriteResultFile("PT", idstring, path, outRGBImgPT);

	m_previousAnalysisOutputDir = path;
	m_idstring = idstring;

	delete progressDlg;

}

GelCorrelationAnalysis::RGBImageType::Pointer GelCorrelationAnalysis::FillRGBImage(GelImage::ImageType::RegionType region, 
											const GelImage::ImageType::SpacingType spacing,
											const GelImage::ImageType::PointType inputOrigin,
											DataImageType::RegionType outputRegion,
											DataImageIteratorType outIter)
{
	RGBImageType::Pointer outRGBImg = createRGBImage(region, spacing, inputOrigin);
	RGBImageIteratorType outRGBIter(outRGBImg, outputRegion);
	outRGBIter.GoToBegin();
	outIter.GoToBegin();
	while(!outRGBIter.IsAtEnd())
	{
		outRGBIter.Set(*(ToRGB(outIter.Get())));
		++outRGBIter;
		++outIter;
	}
	return outRGBImg;
}


void GelCorrelationAnalysis::Permute(std::vector<double>* aVector)
{
	//Implement permutation here.
	int size = aVector->size();
	char out[1024];	
/*
	m_logFile->DisplayText("Before permuted:");
	for(int i = 0; i < aVector->size(); i++)
	{
		sprintf(out, "%f", (*aVector)[i]);
		m_logFile->DisplayText(out);
	}
*/
	for (int i = size - 1; i > 0; i--)
	{
		int random = rand();
		double rand01 = (double)random/(double)(RAND_MAX);
		int index = int(rand01*(double)i + 0.5);
		double tmp = (*aVector)[index];
		(*aVector)[index] = (*aVector)[i];
		(*aVector)[i] = tmp;
		//sprintf(out, "%f %f %d %d %d", tmp, rand01, random, index, i);
		//m_logFile->DisplayText(out);
	}
/*
	m_logFile->DisplayText("Permuted:");
	for(int i = 0; i < aVector->size(); i++)
	{
		sprintf(out, "%f", (*aVector)[i]);
		m_logFile->DisplayText(out);
	}
*/
}

void GelCorrelationAnalysis::NormalizeImage( DataImageType::Pointer outImgD, DataImageType::RegionType outputRegion, DataPixelType aMaxD )
{
	DataImageIteratorType istddev(outImgD, outputRegion);
	istddev.GoToBegin();
	while(!istddev.IsAtEnd())
	{
		DataPixelType val = istddev.Get();
		DataPixelType normalizedVal = val/aMaxD;
		istddev.Set(normalizedVal);
		++istddev;
	}
}

void GelCorrelationAnalysis::WriteResultFile(std::string aImageName, std::string aIdString, gel2de::File* aPath, RGBImageType::Pointer aOutImg)
{
	RGBWriterType::Pointer writer = RGBWriterType::New();
	std::string fileName = aIdString + "_" + aImageName + ".tiff";
	gel2de::File absFile((*aPath), fileName);
	writer->SetFileName(absFile.path().c_str());

	//Flip the image due to ITK/VTK flipping...
	RGBFlipImageFilter::Pointer imageFlipper = RGBFlipImageFilter::New();
	itk::FixedArray<bool, 2> axes;
	axes.SetElement(1, TRUE);
	axes.SetElement(0, FALSE);

	imageFlipper->SetFlipAxes(axes);
	imageFlipper->SetInput(aOutImg);


//	writer->SetInput(aOutImg);
	writer->SetInput(imageFlipper->GetOutput());
	writer->Update();
	writer->Write();
}


void GelCorrelationAnalysis::SetColorTransferFunction(vtkColorTransferFunction *ctf) 
{
  m_ctf = ctf;
}

gel2de::File* GelCorrelationAnalysis::GetPreviousAnalysisOutputDir()
{
	return m_previousAnalysisOutputDir;
}

std::string GelCorrelationAnalysis::GetPreviousIdString()
{
	return m_idstring;
}

GelCorrelationAnalysis::DataImageType::Pointer GelCorrelationAnalysis::getCImage()
{
	return m_outImgC;
}

GelCorrelationAnalysis::DataImageType::Pointer GelCorrelationAnalysis::getSImage()
{
	return m_outImgS;
}

GelCorrelationAnalysis::DataImageType::Pointer GelCorrelationAnalysis::getPImage()
{
	return m_outImgP;
}

GelCorrelationAnalysis::DataImageType::Pointer GelCorrelationAnalysis::getPTImage()
{
	return m_outImgPT;
}

GelCorrelationAnalysis::DataImageType::Pointer GelCorrelationAnalysis::getDImage()
{
	return m_outImgD;
}

GelCorrelationAnalysis::DataImageType::Pointer GelCorrelationAnalysis::getCSImage()
{
	return m_outImgCS;
}

GelCorrelationAnalysis::DataImageType::Pointer GelCorrelationAnalysis::getCSDImage()
{
	return m_outImgCSD;
}

//Non-GUI version, uncomment this to compile as console running.

/*
#include <itkTDistribution.h>

int main()
{

	std::cout << " hei" << std::endl;
	SpearmanRankOrderCorrCoeff *c = new SpearmanRankOrderCorrCoeff(100);

	for(int i = -100; i < 100; i++)
	{

		double rs = 0.01*i;
		double N = 7;
		double t = rs * sqrt((N - 2.)/(1. - rs*rs));

		double probrs = c->BetaI(0.5*(N - 2.), 0.5, (N - 2.)/((N - 2.) + t*t));
		
		double ttest;
		if(rs < 0)
			ttest = 2*itk::Statistics::TDistribution::CDF(t, N - 2.);
		else
			ttest = 2* (1. - itk::Statistics::TDistribution::CDF(t, N - 2.));

		std::cout << "r_s : " << rs << " " << t << " : NR - ITK : " << probrs << " - " << ttest << " = " << probrs - ttest << std::endl;
	}
*/
/*
	GelCorrelationAnalysis *analysis = new GelCorrelationAnalysis();
	analysis->FromXML("gel.xml");
	analysis->RunAnalysis();
	*/
//}



